clear;
clc;
warning off
%% 参数初始化
narvs = 7; % 变量个数
T0 = 20;   % 初始温度
T = T0; 
maxgen = 200;  % 最大迭代次数
Lk = 100;  % 每个温度下的迭代次数
alfa = 0.9;  % 温度衰减系数
x_lb = [0,0,0,0,0,0,0]; % x的下界
x_ub = [20,20,20,20,20,20,20]; % x的上界
% CRA,CL,CLA,C0,C1,CA

%%心室弹性函数信号
% tc = 0.8;
% Tmax = 0.2+0.15*tc;
% t = 0.005:0.005:0.8;
% tn = t/Tmax;
% Entn = 1.55*((tn/0.7).^1.9./(1+(tn/0.7).^1.9))./(1+(tn/1.17).^21.9);
% Emaxleft = 6.0;
% Eminleft = 0.0082;
% Emaxright = 0.00042;
% Eminright = 0.00003;
% Eleft_t= (Emaxleft-Eminleft)*Entn + Eminleft;
% Eright_t = (Emaxright-Eminright)*Entn + Eminright;
% CLV_t = 1./Eleft_t;
% CRV_t = 1./Eright_t;
%% 读取临床数据
load("E:\CS\EECP\数据\a\aRest.mat","Pao");
target_Pao = mean(Pao);
Pao_agj = downsample(reshape(repmat(Pao,1,10).',[],1),16);
%%
Pao_agj = reshape(repmat(Pao_agj,1,20),[],1);
Pao_agj = [Pao_agj;zeros(200,1)];
%%  随机生成一个初始解
x0 = zeros(1,narvs);
for i = 1: narvs
    x0(i) = x_lb(i) + (x_ub(i)-x_lb(i))*rand(1);    
end
%%
load_system('PLM_optimizer1')
set_param('PLM_optimizer1/人体血液循环系统集中参数模型/CRA','Capacitance',num2str(x0(1)))
%get_param('PLM_optimizer/人体血液循环系统集中参数模型/CRA','Capacitance')%0.5945
set_param('PLM_optimizer1/人体血液循环系统集中参数模型/CL','Capacitance',num2str(x0(2)))
set_param('PLM_optimizer1/人体血液循环系统集中参数模型/CLA','Capacitance',num2str(x0(3)))
set_param('PLM_optimizer1/人体血液循环系统集中参数模型/C0','Capacitance',num2str(x0(4)))
set_param('PLM_optimizer1/人体血液循环系统集中参数模型/C1','Capacitance',num2str(x0(5)))
set_param('PLM_optimizer1/人体血液循环系统集中参数模型/CA','Capacitance',num2str(x0(6)))
set_param('PLM_optimizer1/ErightGain','Gain',num2str(x0(7)))
sim('PLM_optimizer1',[0,20]);

sim_Pao_data = simout1.Data();
sim_Pao_data = [sim_Pao_data;zeros(200,1)];
%%
y0 = ClacTargetValue(sim_Pao_data,Pao_agj,1)

%% 定义一些保存中间过程的量，方便输出结果和画图
best_x = x0;
best_y = y0;     % 初始化找到的最佳的解对应的函数值为y0
BESTY = zeros(maxgen,1); % 记录每一次外层循环结束后找到的best_y (方便画图）
%% 模拟退火过程
for iter = 1 : maxgen  % 外循环, 我这里采用的是指定最大迭代次数
    for i = 1 : Lk  % 内循环，在每个温度下开始迭代
        fprintf("外循环：%d  内循环：%d\n",iter,i);
        % 产生随机新解
        y = randn(1,narvs);  
        z = y / sqrt(sum(y.^2)); 
        %z = y;
        x_new = x0 + z*T; 
        
        %对新解的位置进行调整
        for j = 1: narvs
            if x_new(j) < x_lb(j)
                r = rand(1);
                x_new(j) = r*x_lb(j)+(1-r)*x0(j);
            elseif x_new(j) > x_ub(j)
                r = rand(1);
                x_new(j) = r*x_ub(j)+(1-r)*x0(j);
            end
        end
        x1 = x_new;    
        set_param('PLM_optimizer1/人体血液循环系统集中参数模型/CRA','Capacitance',num2str(x1(1)))
        %get_param('PLM_optimizer/人体血液循环系统集中参数模型/CRA','Capacitance')%0.5945
        set_param('PLM_optimizer1/人体血液循环系统集中参数模型/CL','Capacitance',num2str(x1(2)))
        set_param('PLM_optimizer1/人体血液循环系统集中参数模型/CLA','Capacitance',num2str(x1(3)))
        set_param('PLM_optimizer1/人体血液循环系统集中参数模型/C0','Capacitance',num2str(x1(4)))
        set_param('PLM_optimizer1/人体血液循环系统集中参数模型/C1','Capacitance',num2str(x1(5)))
        set_param('PLM_optimizer1/人体血液循环系统集中参数模型/CA','Capacitance',num2str(x1(6)))
        set_param('PLM_optimizer1/ErightGain','Gain',num2str(x1(7)))
        sim('PLM_optimizer1',[0,20]);

        sim_Pao_data = simout1.Data();
        % s2 = sim_Pao_data;
        % s1 = Pao_agj
        y1 = ClacTargetValue(sim_Pao_data,Pao_agj);
       
        if y1 < y0    % 如果新解函数值xiao于当前解的函数值
            x0 = x1; 
            y0 = y1;
        else
            % 根据Metropolis准则判断是否接受这个解
            p = exp(-(y0 - y1)/T); 
            if rand(1) < p  
                x0 = x1; 
                y0 = y1;
            end
        end
        % 判断是否有更好值,更新解


        if y0 < best_y  
            best_y = y0; 
            best_x = x0;
        end
        fprintf("当前MSE = %f,最优MSE = %f\n",y1,best_y);
        disp("当前参数:")
        disp(x1)
    end
    BESTY(iter) = best_y; 
    T = alfa*T;   
    pause(0.01)  
end
%% 输出
disp('最佳的位置是：'); disp(best_x)
disp('此时最优值是：'); disp(best_y)
figure;
plot(BESTY)
%% 设置测试参数
% x0 = [0,0,0];
% f = @fun2
% x = simulannealbnd(f,x0,x_lb,x_ub)
% disp('最优值：');disp(f(x))

set_param('PLM_optimizer1/人体血液循环系统集中参数模型/CRA','Capacitance',num2str(best_x(1)))
%get_param('PLM_optimizer/人体血液循环系统集中参数模型/CRA','Capacitance')%0.5945
set_param('PLM_optimizer1/人体血液循环系统集中参数模型/CL','Capacitance',num2str(best_x(2)))
set_param('PLM_optimizer1/人体血液循环系统集中参数模型/CLA','Capacitance',num2str(best_x(3)))
set_param('PLM_optimizer1/人体血液循环系统集中参数模型/C0','Capacitance',num2str(best_x(4)))
set_param('PLM_optimizer1/人体血液循环系统集中参数模型/C1','Capacitance',num2str(best_x(5)))
set_param('PLM_optimizer1/人体血液循环系统集中参数模型/CA','Capacitance',num2str(best_x(6)))
set_param('PLM_optimizer1/ErightGain','Gain',num2str(best_x(7)))
disp("设置完成")
%% 测试
sim('PLM_optimizer1',[0,20]);
sim_Pao_data = simout1.Data();
%%
% s2 = sim_Pao_2T;
% s1 = Pao_agj
ClacTargetValue(sim_Pao_data,Pao_agj,1)

%% 目标函数
function y = fun(x)
    y = 5*sin(x(1)) - 6*cos(5*x(2)) + 3*x(3);
end
function y = fun2(x)
    y0 = 5*sin(x(1)) - 6*cos(5*x(2)) + 3*x(3);
    y = -y0;
end
function error = ClacTargetValue(s2,s1,display)
    if (nargin<3)
        display = 0;
    end
    
    [C21,lag21]= xcorr(s2,s1);
    C21 = C21/max(C21);
    [~,I21] = max(C21);
    t21 = lag21(I21);
    while t21<0
        t21 =t21 + 100;
    end
    s2_adj = s2(t21+1:end) ;
    s1_adj = s1(1:length(s2_adj))
    mse_t = mse(s2_adj-s1_adj);
    maxvalue_error = abs(max(s2_adj)-max(s1_adj));
    error = 0.5*mse_t + 0.5*maxvalue_error;
    if display
        figure;
        plot(s2_adj,'Color',[1 0 0]);%红色
        hold on;
        plot(s1_adj,'Color',[0 0 1]);%绿色 
    end
    
end
